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Abstract 

Background: Human embryonic stem cell derived cardiomyocytes (hESC-CMs) hold 
high potential for basic and applied cardiovascular research. The development of a 
reliable simulation platform able to mimic the functional properties of hESC-CMs 
would be of considerable value to perform preliminary test complementing in vitro 
experimentations. 

Methods: We developed the first computational model of hESC-CM action potential 
by integrating our original electrophysiological recordings of transient-outward, 
funny, and sodium-calcium exchanger currents and data derived from literature on 
sodium, calcium and potassium currents in hESC-CMs. 

Results: The model is able to reproduce basal electrophysiological properties of 
hESC-CMs at 15-40 days of differentiation (Early stage). Moreover, the model 
reproduces the modifications occurring through the transition from Early to Late 
developmental stage (50-1 10, days of differentiation). After simulated blockade of 
ionic channels and pumps of the sarcoplasmic reticulum, Ca 2+ transient amplitude 
was decreased by 12% and 33% in Early and Late stage, respectively, suggesting a 
growing contribution of a functional reticulum during maturation. Finally, as a proof 
of concept, we tested the effects induced by prototypical channel blockers, namely 
E4031 and nickel, and their qualitative reproduction by the model. 

Conclusions: This study provides a novel modelling tool that may serve useful to 
investigate physiological properties of hESC-CMs. 

Keywords: Embryonic stem cells, Computer simulation, Action potential, 
Pharmacology 



Background 

Human Embryonic Stem Cells (hESCs) are pluripotent cells derived from the blasto- 
cyst stadium of human embryos, having the potential to differentiate in all the three 
embryonic germ layers [1]. Many studies have been carried out to identify the most ad- 
vantageous strategies to drive the differentiation towards the desired cell phenotypes 
thus allowing valuable investigations in basic research and suggesting useful perspec- 
tives for regenerative purposes. In the cardiovascular field, hESCs provide a powerful 
tool to clarify key developmental steps of cardiac embryogenesis [2], to develop reliable 
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in vitro models for drug toxicity screening [3] and are considered a promising source 
for cell-based therapies in pathologies such as myocardial infarction or pace-maker 
center dysfunction [4]. Studies involving hESCs differentiated toward the cardiac 
phenotype are rather demanding due to difficulties such as i) low efficiency process of 
differentiation [5]; ii) dishomogeneity of cell phenotypes; iii) laborious phenotypic 
characterization, e.g. via patch-clamp or multicellular and multi-electrode array record- 
ings [6-8]. Another complication arises from the observation that hESC-derived cardio- 
myocytes (hESC-CMs), like fetal cardiomyocytes (CMs), are electrophysiologically 
immature [6]; their properties evolve during in vitro culturing [6], a phenomenon 
which appears to be regulated by interactions with non-cardiomyocytes in embryoid 
bodies (EBs) [9]. 

To date extensive information is available mostly as unsystematic mass of basic elec- 
trophysiological properties of different hESC lines differentiated toward the cardiac 
lineage and on the modifications occurring during maturation or upon exposure to dif- 
ferent drugs or chemicals. Nonetheless, no attempt was done to systematize current 
knowledge to fully evaluate the impact of individual key ionic currents and of 
excitation-contraction coupling mechanisms on basic physiology in hESC-CMs [10-13]. 

Computational modelling represents a consolidated approach in cardiac research to 
simulate the electrophysiology of single cell or cell-made tissue [14] and the modifica- 
tions induced by chemicals and drugs. This approach usually complements in vitro and 
in vivo experimentation to create a compelling tool able to predict physiological 
responses, abnormal reactions to drug application and to formulate new hypotheses. 

The aim of this work is to develop a computational model of action potential (AP) of 
Hl-hESC-CMs allowing to follow the maturation process. Due to the shortage of mea- 
surements, especially at advanced developmental stages, this model can help to infer 
developmental mechanisms not obvious from the bare measurements. 

Data on membrane ionic currents for this cell line coming from our original mea- 
surements of transient outward potassium, funny, sodium-calcium exchanger currents, 
and data from literature on sodium, calcium and potassium currents in hESC-CMs 
were integrated into the model. To take into account the presence of non-cardiac cells 
in intact EBs, a further model assessment is proposed by coupling the hESC-CM model 
with modelled fibroblasts and evaluating their impact on the AP. The formulated 
model simulates i) the main basic AP features and ii) the developmental changes docu- 
mented during in vitro maturation. 

Methods 

Methods for hESC-CM culture, differentiation and electrophysiological recording are 
described in Additional file 1. hESC-CM were included in the Early or Late group 
according to differentiation time, i.e. from 15-40 and 50-110 days, respectively [6]. 

Model of hESC-CM AP and its transition from Early to Late stage of development 

The starting point in developing the hESC-CM model was a modified version [15] of 
the TenTusscher model of human adult ventricular CM [16], this parent model was 
then largely modified by changing the formulation of almost all the currents to incorp- 
orate all the available data on hESC-CMs and by adding two currents (If and Icar) that 
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are not present in the adult ventricle. Following the classical Hodgkin-Huxley formula- 
tion [17], the cell electrophysiological behaviour is described by Eq. 1: 

dV Ii on 

where V is the membrane potential, C m the membrane capacitance and l ion the sum of 
all the membrane currents. Details on each current are in the following subsections. 

Properties of ion currents based on our recordings or derived from literature data on 
hESC-CMs were integrated into the model to reproduce Early and Late hESC-CM APs. 
Where data from hESC-CMs were not available, observations in ESC-derived or embry- 
onic CMs from different species were considered. Although ionic channels undergo 
complex regulation at a transcript level, the I/V relationship of most currents does not 
change among different developmental stages [18-22]. Hence, we assumed that devel- 
opmental changes in each current, I^, are determined mainly by its quantitative 
change, which can be represented by setting a variable fraction (ratio, Ral^) of the 
current maximal conductance in the adult model. Table 1 summarizes the maximal 
conductance values for the main currents in the model. 

Sodium current (l Na j 

Our I Na formulation slightly changes the original adult model (Additional file 1: Eqs. 
S3, S6 and S8, in Additional file 1). The steady-state inactivation was changed accord- 
ing to Satin et al. data on inactivation dynamics of H9.2-hESC-CMs at the Early stage 
of differentiation (Figure 1A) [23]. 



Table 1 Main developmental changes of ionic currents 





Parameter (units) 


Early 


Late 


Adult [16] 


Species [Ref.] 


l,o 


G to (S/F) max conductance 


1 9.2929 


48.702 


294 


human, EXP, t; human, MOD, [16] 


If* 


Gkt (S/F) max conductance 
Gf (S/F) max conductance 


288.0 
49 


134.4 
20.913 


96 


rat, guinea pig, EXP, [24]; human, 
MOD, [16] 

human, EXP, [6] 


Iki 


G K1 (S/F) max conductance 


240.523 


1154.508 


5405 


human, EXP, [6]; human, MOD, [16] 


IcaL 


G CaL (dm 3 /(F-s)) permeability 


0.0438 


0.0739 


0.175 


mouse, EXP, [24]; human, EXP. 
[6]human, MOD, [16] 


IcaT 


G CaT (S/F) max conductance 


45.8 


9.16 




human, EXP, t; rabbit, MOD, [27] 


iNaCa 


ImaxNaCa (A/F) max current 


17500 


18240 


1000 


human, EXP, t; human, MOD, [16] 


^a 


G Na (S/F) max conductance 


563.844 


14838 


14838 


mouse, EXP, [24]; human, MOD, [16] 


Iks 


Gfe (S/F) max conductance 


15./ 


15.7 


157 


human, EXP, [32]; human, MOD, [16] 


'NaK 


ImaxNaK (A/F) max current 


0.9534 


1.1305 


1.362 


mouse, EXP, [34]; human, MOD, [16] 


IpCa 


Gpca (S/F) max conductance 


0.825, § 


0.825, § 


25 


human, MOD, [16] 


'up 


Uaxup (mM/s) max flux 


0.0565 


0.1403 


0.425 


mouse, EXP, [24]; human, MOD, [16] 


Irel 


UaxRei (mM/s) max flux 


0.274 


9.88 


24.7 


mouse, EXP, [24]; human, MOD, [16] 


'leak 


LaxLeak (1/s) max rate 


0.0004 


0.024 


0.08 


mouse, EXP, [24]; human, MOD, [16] 


IbCa 


G bCa (S/F) max conductance 


0.118, § 


0.592, § 


0.592 


human, MOD, [16] 



Maximum conductances, currents and fluxes for Early and Late human embryonic stem cell-derived and adult ventricular 
cardiomyocyte models. The species considered for the Early and Late date are reported. Adult data refers to the 
TenTusscher model [16]. t new measurements and data; § values chosen to reproduce at best the AP shape; * not only 
the maximal conductance/current was changed but also the current formulation (see Additional file 1). EXP: experimental 
data; MOD: used for modelling. 
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Our model 
Exp. data [23] 
Ten Tusscher 
model [16] 




A"? A? & & 
* C cfcfcfcfcf 
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Cav 3.2 

a -acti ni n 
GAPDH 




5^ 

Q. 

5z -6 



-Our model Early 
Early exp. 
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-12-1 

-150 



Early 



Late 



-100 
V(mV) 



Q. 

5 




Early 



-50 

V(mV) 



Figure 1 Currents: experimental data and model fitting. (A) l Na inactivation for Early human embryonic 
stem cell-derived cardiomyocytes (hESC-CMs). (B) Ic 3 l normalized current-voltage (l/V) curve (Late stage). 
(C) Cav3.1, Cav3.2 and GADPH expression (l CaT ) for HI -hESC-CMs. (D) l to l/V curves. (E) l f l/V curves and HCN 
quantitative expression {inset). (F) l NaCa elicited by a voltage ramp protocol 



As also reported in [24] for rodent ESC-CMs we considered a small expression 
of I Na at the Early stage and full expression at the Late one; the expression at the 
Early stage was further reduced with respect to [24] (from 0.08 to 0.038) in order 
to fit properly our experimental maximal upstroke velocity (Vmax) and the action 
potential duration (APD). 

RaI Na = 0.038 (Early); Ral Na = 1 (Late) 
L-type calcium current (l Ca J 

Due to lack of specific data on I CaL at the Early stage, we tuned the permeability, 
the Ca 2+ dependent inactivation gate, f Ca (Additional file 1: Eqs. S25-S31), the time 
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constant of the voltage dependent inactivation gate, Tf (Additional file 1: Eq. S33), 
and the steady-state activation gate, d^ (Additional file 1: Eq. S19) to reproduce 
the experimental features of AP, in particular APD. 

We then slightly modified d m in order to fit our experimental recordings of Ic a L 
at Late developmental stage (Figure IB, Late stage, n = 1), whereas inactivation 
parameters became equal to those of adult CMs in the Late stage formulation. We 
maintained the ratio between Late and Early conductances proposed in [24] based 
on data in mice and guinea pigs: 

Ralcai = 0.25 (Early); RaI Ca L = 0.422 (Late) 

T-type calcium current (l CaT j 

At variance with the adult model, we included I CaT , on the basis of different ex- 
perimental evidence. First, I CaT was reported to be highly expressed during fetal 
heart development and gradually decline after birth, becoming restricted to the 
conduction and pacemaker cells [25]. Secondly, I CaT is functionally expressed in 
mouse ESC and is downregulated during cardiac differentiation: I CaT channel subu- 
nits Cav3.1 and Cav3.2 expression decreased to approximately 46% and 24%, re- 
spectively, at 23.5 days of differentiation with respect to 9.5 days [26]. Finally, our 
qualitative RT-PCR measurements in Hl-hESC-CMs show a clear expression of 
Cav3.1 and Cav3.2 at both stages (Figure 1C). We used the I CaT formulation pro- 
posed in [27] in their sinoatrial node cell model, with progressively decreasing 
scaling factors for the maximal conductance: 

RaI CaT = 0.25 (Early); RaI CaT = 0.05 (Late) 

Transient outward (l to ), rapid and slow delayed rectifier (l Kr , l Ks ) IC currents 
I t o properties were based on original data obtained in Hl-hESC-CMs. Figure ID 
shows the I/V relationship for peak K + currents evoked by depolarizing steps in 
Early and Late CMs (median data and interquartile, n = 10 Early stage and n = 9 
Late stage). A 10 mV positive shift was applied to experimental data to account 
for the use of Cd 2+ to block I CaL , as done by [28]. According to our previous 
observations [6], I to activation properties are similar to those described in native 
cardiac cells. Maximum conductances and steady-state activation were calculated 
by fitting experimental data (Figure ID and Additional file 1: Eq. S43): 

Ral t0 = 0.065622 (Early); Ral t0 = 0.165653 (Late) 

In accordance with various in vivo and in vitro experimental data in fetal 
guinea pigs [20] rats [29,30], and mice [31], I Kr maximal conductance was greater 
than in the adult CMs and it decreased during maturation, we chose the specific 
expression ratios in order to mimick maximum diastolic potential (MDP) and the 
APD at the different developmental stages: 

RaI Kr = 3 (Early); 7? aI Kr =1.4 (Late) 

I Ks conductance in the Early developmental stage was set in order to achieve a 
good fitting of the I/V curve (Additional file 1: Figure SI) obtained by [32] on 
Early hESC-CMs. On the basis of data on embryonic murine heart [19,24], where 
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Early and Late developmental stages seem to share the same Ik s conductance, a 
single value of RaI Ks was used: 

Ralxs = 0.1 (Early and Late); 
Inward rectifier IC current (l K1 ) 

In hESC-CMs I K1 is very small but not absent at the Early stage, then it increases 
during the development, as also reported for rat and guinea pig [24]. Conductances 
were identified according to our previous data [6], showing a ratio between Late 
and Early current density at -90 mV of 5.42. In order to reduce the MDP we set 
in our model a ratio of 4.8, within the variability of our data and introduced a 
shift of the voltage dependence for the inward rectification factor, x kloo (Additional 
file 1: Eqs. S67-S69), without modifying the current reversal potential. 

RaI KX = 0.0445 (Early); RaI K1 = 0.2136 (Late) 
Funny current (l f ) 

A key step of the formulation of a specific hESC-CM model consisted in integrat- 
ing the hyperpolarization-activated cyclic nucleotide-gated or funny current, that is 
reported to be one of the main contributor for the spontaneous beating of pace- 
maker cells [33] and hESC-CMs [6]. I f formulation (Additional file 1: Eqs. S71-S73) 
and conductance were obtained by fitting recordings performed on Early Hl- 
hESC-CMs (Figure IE, Early stage, n = 4). 

For the Late stage, in accordance to our previous data [6], we assumed that current 
density decreased over maturation to an extent equal to the drop of cumulative HCN 
transcript expression (Figure IE, inner panel): the estimated ratio between Early and 
Late stage was 2.34. 

Ralf = 0.5389 (Early); Ral f = 0.23 (Late) 
Sodium-potassium pump (l NaK ) and sodium-calcium exchanger (l Na c a ) 

Since data on hESC-CM I NaK were not available, maximal current density was set tak- 
ing into account its influence on the diastolic depolarization rate (DDR) and frequency 
of spontaneous beating (F) and reflecting the maturation related growth of INaK ex- 
pression according to experiments by [34] on mouse ESC-CMs: 

RaI NaK = 0.7 (Early); RaI NaK = 0.83 (Late) 

As far as I^Ca is concerned, we used original experimental data from Hl-hESC-CMs. 
Voltage ramp (from -120 to + 50 mV) protocols elicited an almost linear ^aCa I/V rela- 
tionship (Figure IF, n = 6) that showed an inward and outward mode at both develop- 
mental stages. Fitting of experimental data led to modify the original maximal current 
density and the extra factor a in the INaCa expression [16]: 

Ralmca = 17.50 (Early); RaI NaCa = 18.24 (Late) 
a = 0.8 (Early), a = 0.38 (Late) 

Sarcoplasmic reticulum (SRj currents 

The maximal values for the uptake (I up ), release (I re i) and leakage (Ii ea k) currents were 
tuned to simulate the ryanodine induced reduction of Ca 2+ transient amplitude 
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reported in [11] at the Early and in [12] at the Late stage. Increases in maximal current 
densities at the Late stage were based on the rodent ESC-CM model [24]. 

RaI Up = 0.133 (Early); RaI Up = 0.33 (Late) 
RaI Re i = 0.0111 (Early); RaI Rd = 0.4 (Late) 
RaI Leak = 0.005556 (Early); RaI Leak = 0.3 (Late) 

Sarcolemmal calcium pump, i pCal and background current, l bCa 

Since data about these currents are not available, we chose the following values to re- 
produce at best the AP shape: 

Ralcap = 0.033 (Early and Late); 
Rahca = 0.2 (Early) ; Rahca = 1 (Late); 

Cell capacitance and dimensions 

Median values of measured cell membrane capacitance were used to set cell dimen- 
sions (see Additional file 1). 

Sensitivity analysis 

A sensitivity analysis was performed according to the procedure reported in [35,36], 
opportunely adapted to our model. The main differences with respect to [35] are: 
(i) our hESC-CM model is not stimulated by an external source and (ii) we concentrate 
our analysis on the impact of the ratios Ral xx on the AP shape. One ratio was varied at 
time by -20%, -10%, +10% and +20% respectively. 
Considering the following ratios (parametrs) "p" 

p = {RaI Na , RalcaL, Ralf, RaI to , RaI K i, RaI Kr , RaI Ks , RaI N aK, RaI N aCa} 

and the AP features (characteristics) "c" 

c = {MDP, Vmax, APD30, APD50, APD70, APD90, DDR, F}, 

computed after 300 seconds of simulation (assuming the steady state condition) the in- 
dexes percentage of change (D c _ Pia ), sensitivities (S Cj p j+2 o% and S CjPi _ 20 %) and relative sen- 
sitivities (r Ci p j+2 o% and r CjP> . 2 o%) were calculated as follows: 



\f^p^a ~t~ ^control ) 

c,p,c ~ 



D c ,p,a = L 100 

^control 



D c ,p,+20% „ D c ^_20% 



J c,p,+20% — q 2 ' °c,/> -20% q 2 

Scj>,+20% S cp _2Q% 

r c,p,+20% — Yc 1 ' r c,p -20% — 1 

\\,p,+20% | maXjC \\p,-20% I max _ c 

Splitting the original S c p and r c p [35] was necessary since several tests resulted in no 
spontaneous APs, thus making impossible calculating the AP features and all their D c pa . 
However, for each ratio at least one D CjP _ a (D CiP +2 o% or D CiPi . 20 %) was available thus allow- 
ing to get the asymmetrical sensitivities. 

Interaction with in silico fibroblasts 

In order to preserve intracellular milieu and cell-to-cell communication, AP recordings 
were not performed on single cells but on EBs, aggregates containing different cell 
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phenotypes among which hESC-CMs and fibroblasts. To test the interaction between 
these kinds of cells and assess the effect on AP simulation, an additional mammalian 
fibroblast model, resistively coupled to the hESC-CM, was developed according to [37- 
39]. To this aim, the hESC-CM membrane potential equation was modified as follows: 



where V flbro is the fibroblast potential, G gap is the conductance of the hESC-CM-fibroblast 
coupling (G gap = 1 nS) and N f is the number of coupled fibroblasts. Details on the fibro- 
blast model are reported in the Additional file 1 (Additional file 1: Eqs. S92-S94). 

Drug simulations 

To reproduce the inhibition of SR Ca 2+ release and the consequent Ca 2+ transient re- 
duction in ryanodine experiments the simulation was performed zeroing I up and I re i. 
To simulate E4031 (I Kr blocker) and nickel (Ic a T» Ii<r and iNaCa blocker) effects, in steady 
state conditions, a step reduction of conductance for the targeted currents was imple- 
mented in the model. The amount of conductance reduction for each current, which is 
reported in the Results Section, was chosen within the range of blocking action of the 
drug in order to better reproduce the specific experimental result obtained on a single 
cluster of hESC-CMs. 

Numerical implementation 

Differential equations were implemented in MATLAB (The Math Works, Natick, MA) 
and solved using odel5s. 

Results 

hESC-CM model 

Figure 2A shows simulated AP profiles for Early hESC-CMs obtained using our model: the 
modifications introduced with respect to the adult ventricular model were sufficient to elicit 
spontaneous beating. A comparison with experimental Early AP profiles obtained on intact 
EBs by multicellar recordings is provided in Table 2: a global comparison was done by cal- 
culating typical morphological parameters (AP features) for both experimental and simu- 
lated APs. This analysis demonstrated that our hESC-CM model was able to reproduce 
most of the experimental AP features, including AP duration (APD) at 30 (APD 30 ), 50 
(APD 50 ), 70 (APD 70 ) and 90% (APD 90 ) of repolarization, Vmax, F and the DDR. Simulated 
and experimental data differed for the MDP, which is more negative in the simulation. 

Transition from Early to Late stage of development: effects of maturation 

The simulated AP profile at the Late stage and the comparison with experimental APs 
obtained on intact EBs by multicellar recordings are reported in Figure 2B and Table 2, 
respectively. These results show that the changes introduced in the model parameters 
between the Early and Late stage allow to reproduce the documented [6] maturation 
effects on AP shape. In particular, during the transition from the Early to the Late 



dV 
dt 



lion H - 'Igap 



in 



(2) 




(3) 
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20 

> o. 
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Time (s) 
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234567 89 10 
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Figure 2 Action Potential simulations. Simulated action potentials for Early (A) and Late (B) hESC-CMs clusters 



stage, APD and Vmax increase while spontaneous rate and slope of diastolic 
depolarization decrease. 

Spontaneous firing and action potential shape dependence on current conductances 

The sensitivity analysis performed on the Early and Late models was aimed to assess 
how variations in the maximum conductances of the most important membrane cur- 
rents affect (i) the phenomenon of spontaneous beating and (ii) the AP shape. 

The spontaneous firing activity was not triggered at the Early stage in 2 tests only: 
RalcaL -20% and RaI NaCa +20%. At the Late stage the spontaneous activity showed to 



Table 2 Action potential features of experimental and simulated action potentials 

Early Late Delta % 





Exp 


Sim 


Exp 


Sim 


Exp 


Sim 


APD 30 (ms) 


117 (1074132) 


121 


163 (844256) 


165 


+ 40 


+ 36 


APD 50 (ms) 


169 (149-=- 184) 


163 


258 (1484374) 


224 


+ 53 


+ 37 


APD 70 (ms) 


205 (1814-219) 


191 


325 (1794465) 


283 


+ 59 


+ 48 


APD 9 o (ms) 


229 (2164-241) 


225 


419 (2234506) 


350 


+ 83 


+ 56 


Vmax (mV/s) 


4610 (330845612) 


4123 


5566 (519847047) 


5620 


+ 42 


+ 36 


MDP (mV) 


-55 (-624-37) 


-76 


-51 (-624-42) 


-73 


-7 


- 4 


DDR (mV/s) 


14.4 (11.14284) 


9.7 


9.9 (8.0411.1) 


7.1 


- 31 


- 27 


F (bpm) 


27 (23453) 


29 


22 (20429) 


24 


- 17 


- 17 



Delta %, Late VS Early; APD, action potential duration; V max maximum upstroke velocity; MDP maximum diastolic 
potential; DDR diastolic depolarization rate; F frequency. Data reported as median and interquartile. All the experimental 
action potential features were measured on hESC-CM AP at the Early and Late stage [50]. 
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be sensitive to more currents: the firing activity was triggered but stop after few dozens 
of seconds of simulations for RaI Na +20%, RaI Ca L-20%, Ral f -20%, RaI K1 +20%, RaI Kr 
+20%, RaI NaK +20% and RaI NaCa +20%. 

When the spontaneous beating allowed to reach a steady state condition, the AP fea- 
tures absolute/relative sensitivities were computed (Figure 3) as reported in the Meth- 
ods section. At the Early stage MDP is affected by the inward IcaLi iNaCa (inward 
during the late-repolarization) and the outward I Kr - The effect of If is extremely small, 
since at the MDP potential it is not activated yet, due to its high time constant. Also 
the I K i effect is small since at this developmental stage the I K1 expression is low. The 
Late stage shows a more stable MDP (maximum variations: 7-8%), but the I K1 effect is 
stronger, as consequence of a maturation-related increment in I K1 expression. As 
expected, Vmax is mainly affected at both stages by the inward currents mainly acting 
during the upstroke: I Na and I CaL , while the rate of spontaneous beating resulted to be 
more sensitive to I CaL , I NaK , I N aCa> If (Early stage)and I K1 (Late stage). I CaL , I NaK and I K1 
had the most strong effect on DDR and F at both stages; these AP features show also 
an important sensitivity to If increments at the Early stage. It is interesting to note that 
also iNaCa reduction increase DDR and F at both stages but particularly at the Early 
one: a smaller INaCa (inward during late-repolarization) makes MDP more negative 
(-81 mV vs -76 mV, RaI NaCa -20% and control respectively) allowing a greater activa- 
tion of If and thus increasing both DDR and F. Variations of the outward currents I K i 
and I NaK show the role of these currents in stabilizing the diastolic potential and, at the 
Late stage, increments block the spontaneous beating. APD is mainly affected by the in- 
ward currents Ic a L and INa. by INaCa (outward during the upstroke), INaK (Late only) 
and the outward Ik„ especially during the late-repolarization. At the Early stage, the 
APD decreases after I CaL increments: this counterintuitive result can be explained by 
using the model. In fact, it is due to a higher AP peak (20 mV vs 8.3 mV, RaIc a L + 20% 
vs control): the higher reached membrane potential allows a stronger I Kr activation thus 
reducing APD. At the Late stage this phenomenon is not present, since in the AP peak 
increment is smaller (34.2 mV vs 30.2 mV, RaI CaL + 20% vs control), and I Kr itself is 
smaller due to maturation. This different contribution of inward/outward currents to 
APD likely underlies the diverse APD rate-dependence (Additional file 1: Figure S2). 

Coupling with fibroblasts 

In order to partially overcome the discrepancy in MDP between the model and the ex- 
perimental measurements, we introduced the contribution of fibroblasts, which are an 
essential component of EBs also for CM maturation. To assess the relevance of this 
issue in our model, we considered a simple system composed of a single hESC-CM 
coupled with 1 and 2 fibroblasts for each stage. Changes in basal AP due to fibroblast 
coupling are reported in Figure 4A and 4C, for the Early and Late phase, respectively, 
while a quantification of these changes is summarized in Figure 4B and 4D. Changes in 
most of the AP features as a function of fibroblast number were similar between the 
two stages. In particular, we observed an increment in DDR and rate while AP ampli- 
tude (APA) decreased. The effect on APD was different in the two stages. In the Early 
hESC-CM, lacking the plateau phase, the major effect is an inward current flowing 
from the coupled fibroblasts into the hESC-CM during the late repolarization phase, 
when hESC-CM membrane potential is negative to the fibroblast resting potential, 
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Figure 3 Relative sensitivity maps for the Early (A) and Late (B) models. AP features (rows) VS ratios 
Ral xx (columns) used to rescale the maximum conductances/currents/fluxes were considered. For each Ral^ 
the relative sensibilities at -20% and + 20% were taken into account. White color indicates maximum 
relative sensitivity of a particular AP feature among all ratios, whereas black indicates AP feature and ratio 
are independent. White X indicates the absence of spontaneous firing during a particular test. Percentages 
in each white box indicate the maximum absolute sensitivity of the AP feature correspondent to that row 
for all ratios. Negative sign indicates that AP feature and ratio vary inversely 




promoting slowing of repolarization and APD increase. In the Late hESC-CM, showing 
significantly longer AP with respect to Early, the effect of an outward current, flowing 
towards coupled fibroblasts at depolarized potential and promoting hESC-CM repolari- 
zation and APD decrease, is also important. The overall result in Late hESC-CM is a 
decrease of APD 30 and APD 50 whereas APD 70 and APD 90 were almost not affected 
by coupling with fibroblasts. Importantly, as the number of coupled fibroblast 
increases from 0 to 2, a small rising of MDP occurred both at Early (-76 to 
-73 mV) and Late stage (-73 to -72 mV). This effect was accompanied by a re- 
duction of Vmax at Early (4123 to 3443 mV/s) and Late stage (5620 to 2554 mV/ 
s). At the same time, the membrane potential of coupled fibroblasts mimics the 
hESC-CM AP (data not shown), oscillating between -74 and 4 mV at the Early 
and between -72 and 26 mV at the Late stage. 



Intracellular calcium 

To assess the relevance of RyR-mediated SR Ca 2+ release in our model, we simu- 
lated cytoplasmic Ca 2+ oscillations in control conditions and after blockade of Ca 2+ 
release. In control conditions, the Early model showed intracellular Ca 2+ diastolic 
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Figure 4 Electrical coupling between human embryonic stem cell-derived cardiomyocyte and 
fibroblasts. (A, C) Effects on action potentials (APs), Early and Late stage respectively, of a variable number 
(N fl ) of fibroblasts. (B, D) Consequent changes of AP features, Early and Late stage respectively: both stages 
show an increment of the rate of spontaneous beating and an APA reduction. APD is slightly prolonged at 
the Early stage (especially APD70 and APD90) while at the Late stage the fibroblast effects are minor, 
showing only a slight reduction of APD30 



and systolic concentrations of 0.026 ujvl and 0.141 ujvl, respectively, with an ampli- 
tude of the transient of 0.115 uM. The mean rate of decay was 0.123 uM/s while 
the maximum upstroke velocity (VmaxCa, upstroke) was 1 uM/s and the maximum 
decay velocity (VmaxCa.decay) was 0.52 uM/s. Control values for the Late model 
were 0.063 uM and 0.506 uM for diastolic and systolic concentrations, 0.443 uM 
and 0.340 uM/s for amplitude and mean rate of decay. VmaxCa,upstroke was 13 
uM/s and VmaxCa,decay was 1.4 uM/s. 

The blockade of the SR channels and pumps reduced the Ca 2+ oscillation amplitude 
by only 12% at the Early stage, whereas the RyR-induced reduction was larger, 33%, at 
the Late stage (Figure 5 A and 5B). These results are similar to those reported experi- 
mentally in Hl-CMs by [11] after 18-24 days and by [12] in cells assessed 30-40 days 
"post-beating", likely corresponding to 50 days of total time of differentiation. A more 
detailed comparison was performed between the data reported in [11] on caffeine- 
sensitive Hl-CMs and our Early stage, also considering the transient elicited in [11] in 
stimulation condition. In order to compare the experimental and simulated VmaxCa, 
upstroke and VmaxCa.decay in control conditions we normalized them using the 
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amplitude of the transient, since experimental data were reported in terms of fluorescence 
while ours as concentration values. Estimating from [11] transient amplitude 0.16 F340/ 
380 and VmaxCa,upstroke 1.5 F340/380/s we got a normalized VmaxCa,upstroke of 9.3 
1/s. Normalizing our Early VmaxCa,upstroke, our model reproduced 8.7 1/s. Comparing 
in the same manner the VmaxCa,decay (estimated experimental VmaxCa,decay 0.3 F340/ 
380/s) we got 1.9 1/s vs 4.5 1/s, respectively experimental [11] and simulated. The last 
comparison regards the reduced VmaxCa.upstroke value after RyR application: in [11] an 
experimental value of 70% was reported, while our model simulated 78%. 

Drug simulations 

Finally, we performed a preliminar test of the qualitative reproduction of the main effects 
on AP of well-established channel blockers. E4031 is known to prolong cardiac AP 
through I Kr selective blockade. We tested the effect of E4031 on a single spontaneously 
beating cluster of hESC-CMs at Early stage and observed a progressive decrease of the 
pacing frequency, a depolarization of MDP and an increase of APD involving all phases of 
repolarization (Figure 6A, left). These effects were fast on our experimental system and 
after 60 s from drug application the cluster showed a complete stop of its spontaneous 
electrical activity and beating [6] . To simulate this effect in our model, G Kr was decreased 
by 50% of its Early value and this resulted in AP prolongation (Figure 6A right). In fact the 
AP lengthening measured at different values of repolarization were 12% (APD 30 ), 18% 
(APD 50 ), 47% (APD 70 ) and 110% (APD 90 ). In the simulation MDP was also depolarized by 
4%, in accordance with a weaker contribution of repolarizing current, and frequency was 
reduced by 27%. In our model a decrease of I Kr current larger than 50% produced a block 
of spontaneous beating, similarly to what observed experimentally [6]. 

In a single cluster of hESC-CMs at the Late stage, application of E4031 produced a fast 
depolarization of MDP and a remarkable increase in spontaneous rate with characteristic 
small amplitude oscillations of membrane potential (Figure 6B, left). At this stage, simula- 
tion was obtained by 60% reduction of G Kr that resulted in AP modifications qualitatively 
similar, even if smaller, to those observed experimentally. In particular, rate was increased 
by 97%, APD 70 by 19%, APD 90 by 65%, MDP was depolarized by 17% (Figure 6B, right). 

At millimolar concentration nickel is well known to largely block iNaCa an d IcaT- Fur- 
thermore, it has also been reported to block at a very high extent I Kr in sinoatrial node 
cells [40]. All these currents are involved in different phases of spontaneous beating 
generation and AP. In the Early stage, application of nickel produced a marked MDP 
depolarization, a sustained APA reduction of and an increase of spontaneous rate, 
which lasted unaltered over time (Figure 6C, left). Conversely, at the Late stage, nickel 
first slowed down and then blocked completely electrical activity (Figure 6D, left). At 
both stages, simulations were performed reducing I CaT and I Kr by 90% and iKaCa by 75% 
of their initial values. The effects on APA, MDP and beating rate detected at the Early 
stage were well reproduced (Figure 6C, right). Blockade of spontaneous activity 
observed at the Late stage was also clearly mimicked in our model (Figure 6D, right). 

Discussion 

The major aim of this work was to develop a mathematical model able to reproduce 
the basal electrophysiological properties of hESC-CMs and the modifications induced 
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Figure 5 Effects of Ryanodine. Simulation of the reduction of calcium transient amplitude under the 
effect of Ryanodine at Early (A) and Late (B) stages 



by in vitro maturation. This approach, has been applied for the first time to hESC-CMs 
and the model represents the first computational tool for studying the fundamental 
physiology of hESC-CMs, in particular those derived from the HI cell line. Moreover, 
introducing the contribution of fibroblasts we approached the simulation of the proper- 
ties of beating intact EBs, which represent the elementary functional unit able to pro- 
mote electrophysiological maturation of hESC-CMs and therefore suitable for 
screening ion channel blockers and assessing the cardiac safely of drugs. 
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Figure 6 Effects of current blocking with E4031 and Nickel. Single experiments [left) and simulations 
[right) on the effect of 10 uM E4031 (l Kr blocker) and of 2 mM Nickel (l NaCa , tar, In blocker) at Early (A and 
C) and Late (B and D) stages 



We reached the goal of reproducing the spontaneous AP profile of hESC-CMs, 
strongly focusing on the modifications occurring during the transition of intact EBs 
from the Early to the Late developmental stage. Specifically, we reproduced key modifi- 
cations documented experimentally on intact EBs by multicellular recordings, including 
the decrease of spontaneous AP frequency in association with a reduced DDR, meaning 
that the automaticity was slowed down during simulated maturation, similarly to ex- 
perimental observations. The maturation-related increase of APD was also satisfactorily 
mimicked by our model, thus reflecting the modifications identified for many ionic cur- 
rents occurring during in vitro maturation. Of note, we chose this specific electro- 
physiological approach on intact EBs due to several advantages, including the 
preservation of tissue architecture that allows the detection of single cell transmem- 
brane voltage resulting also from the contribution of neighbouring cells. The latter, be- 
side non myocytic cells, may be represented by CMs with similar or different 
phenotypes, with a relative composition of different phenotypes depending mainly from 
the developmental stage under evaluation. In fact Early hESC-CMs have a rather homo- 
geneous phenotype (mostly sinus nodal like) that progresses into different phenotypes 
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in the Late phase (atrial, ventricular and sinus nodal). Therefore, our experimental and 
simulation data readily reflect the degree of phenotype dishomogeneity in most of ac- 
tion potential parameters, resulting from the contribution of individual CMs at differ- 
ent stages (throughout the developmental process) and the fate (ventricular, atrial or 
nodal). To further confirm this statement, the median measurements of action poten- 
tial duration at different percentage of repolarization have a lower interquartile range 
in the Early phase compared to the Late (see Table 2). 

Based on the integration into our model of original and literature data on different ionic 
currents, some considerations can be drawn on their role in the hESC-CM development. 

I to is known to play a key role during the early repolarization phase leading to the 
notched shape and plateau phase present in human adult ventricular and atrial cells 
[41]. In our experimental conditions we found that I to density increased from the Early 
to the Late phase. This result is in line with our previous data demonstrating a func- 
tional and molecular up-regulation of I to during hESC-CMs maturation. It further con- 
firms that this channel is a fundamental marker of functional CM maturation among 
mammalians [6]. These modifications were favourably integrated in our model, thus 
demonstrating their positive contribution to its final formulation. However, our mea- 
surements of I to showed substantial current density also for negative potentials. For this 
reason I to half-maximal activation was shifted from 20 to -5 mV (Figure ID) whereas 
the steepness of the model curve was kept slightly greater than the experimental one, 
as already done and discussed in the adult model [16]. 

A different set of original data integrated in our model is related to If. It is an inward 
current involved in the generation of spontaneous electrical activity, identified and 
characterized in our previous study on Hl-hESC-CMs [6]. In the present study we fur- 
ther characterized this current identifying its functional relevance over an extended 
period of time. We calculated maximal conductance for the Early stage and extrapo- 
lated its value at the Late one on the basis of a cumulative reduction of HCN total tran- 
script. Although functional properties of currents do not uniquely depend on the 
amount of channel transcript/s, we hypothesised that I f density is likely to decrease 
during hESC maturation on the basis of different experimental evidences. First, the rate 
of spontaneous beating and the slope of the diastolic depolarization decrease with mat- 
uration (Table 1); secondly, I f downregulation is a well established marker of functional 
maturation of native CMs [42,43]. These observations led us to integrate in our model 
a If contribution declining over maturation time. As observed for I to , this implementa- 
tion had a positive effect on the model mimicking properties. As shown in Figure 3, 
modulation of If mainly affects in the model the rate of spontaneous beating and the 
DDR at both stages but especially at the Early one. The Late stage represented a sce- 
nario characterized by a smaller I f and which evolves towards cells with no automati- 
city, so the 20% reduction of Ral f led to no spontaneous activity. 

Experimental evidence documents the occurrence of developmental changes of I^Ca 
in different animal models, with maximal current density lowering throughout matur- 
ation [24]. In humans, I^Ca expression peaks at mid-gestation, overcoming that of the 
adult heart [44]. In hESC-CMs, [11] reported that INaCa expression in the Early phase is 
higher with respect to that found in human fetal and adult ventricular CMs. This evi- 
dence is in line with our experimental results (see Figure IF), showing larger current 
density both in Ca 2+ outward and inward modes, at the Early stage with respect to the 
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Late one. By contrast, a recent study [45] performed on a similar model led to opposite 
results, i.e. maximal current occurring in the Late stage. The explanation for such a dif- 
ference is not obvious; diverse developmental window and/or cell phenotypes (ventricu- 
lar and atrial vs ventricular) may account for these discrepancies. Model-based analysis 
showed that at both stages iNaCa (Figure 3) modulation affects basically all the AP fea- 
tures (tests with RaI NaCa -20%) and the occurrence of spontaneous beating (tests with 
RaI NaCa +20%). 

Other relevant maturation-related current modifications introduced in our model are 
consistent with previous observations in rodent ESC-CMs. In particular I Kr decrease 
and I Na increase confirmed to be important to achieve AP changes such as AP pro- 
longation and increase in Vmax observed with hESC-CM development. As expected 
I Na had the strongest influence on Vmax while I Kr had a considerable effect on the 
APD, especially in the latest phases of the repolarization (see Figure 3). 

The sensitivity analysis summarized in Figure 3 helps in assessing the maturation- 
related effects on the spontaneous contractile activity of hESC-CMs, in particular 
showing that at the Late stage it is more sensitive to current variations that at the Early 
stage. At the Early stage the reduction of the depolarizing I CaL (necessary for the up- 
stroke) and the increment of I^Ca (inward current during the late-repolarization) 
caused the spontaneous beating stop. In addition to the reduced I CaL and increased 
I Na c a tests, the Late stage showed no spontaneous beating also for increments of the 
main repolarizing currents I Kr > Iki> Iks an d lNaK> thus supporting the hypothesis of a 
weakening of the spontaneous depolarization during maturation. Moreover, also the ef- 
fect of the 20% If reduction at the Late stage is indicative of this mechanism since If 
was already reduced during the transition from the Early to the Late stage. 

We also tested coupling of fibroblasts with the hESC-CM model. Our hypothesis was 
that such approach could improve the quality of simulated APs that are generated by a 
complex system comprising CMs embedded in a cluster of different cells, such as fibro- 
blasts. Indeed, one specific AP feature predicted by our model (MDP) differed substan- 
tially with respect to the experimental data. At both developmental stages, effects of 
coupling were small and consisted of MDP depolarization and increase of DDR and fre- 
quency. These effects can be explained by observing that during diastole fibroblast 
membrane potential is less negative with respect to that of CMs, causing a small de- 
polarizing current flowing into CMs. Although limited, these modifications collectively 
move our simulation towards the experimental measurements, therefore improving the 
mimicking potential of our model. On this basis, it would be interesting to explore the 
possibility to enlarge the cell network by including different cell phenotypes (e.g. endo- 
thelial cells) possibly present in vivo. Globally, our results are in line with similar stud- 
ies [37], where coupling of the TenTusscher model to fibroblasts led to a slightly 
depolarized resting potential, reduced APA and to the electrotonic modulation of the 
fibroblast potential by the coupled CM. 

Recently, increasing attention has been focused on the development and maturation 
of the SR activity in hESC-CMs. Accumulated evidence points to a significant activity of 
Ca 2+ release from SR [11,12], even if less organized and regulated than in mature CMs. 
In fact, while in the Early phase Hl-hESC-CMs express SERCA2a, other proteins such 
as calsequestrin, triadin and junctin are almost absent [46]. Moreover, at 40-50 days of 
differentiation, T-tubules remain undetected on the sarcolemma, suggesting a 
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topological environment different from mature, where Ca + influx through Ca + channels 
in the T-tubule is tightly associated with sarcoplasmic RyR channels [10]. Overexpression 
of calsequestrin in Hl-hESC-CMs failed to induce the growth of T-tubules even though 
SR load and Ca 2+ transient amplitude became more similar to those of mature CMs [47]. 
Simulations of Ca 2+ transients using our model led to results fully consistent with the ex- 
perimental data reported above. In fact, upon blockade of Ca 2+ release from SR we found 
a decreased Ca 2+ transient amplitude by 12% in the Early, and by 33% in the Late phase, 
similarly to the results reported experimentally [11,12]. Overall, these results indicate that 
in the Early stage Ca 2+ cycling is mainly governed by sarcolemmal fluxes, while upon mat- 
uration Ca 2+ release from SR increases its contribution moving toward a functional inte- 
gration with sarcolemma to generate rhythmic activity. This evidence is in line with 
recent data obtained in late-stage mouse ESC [46] and further extends the predictive po- 
tential of our model to intracellular Ca 2+ handling processes. Simulated block of Ca 2+ re- 
lease also increased the oscillation frequency (see Figure 5); we are not aware of 
experimental reports on this specific aspect in hESC-CMs which, on the other hand, 
seems in contrast with the effect of RyR knock out in ESC-CMs [48]. However, in the 
same conditions, a consistent increase in the frequency of AP-induced Ca 2+ transients is 
apparent also in [12] (Figure 6), thus suggesting that further investigation supported also 
by model based analysis might provide mechanistic insights on this issue. 

In our model, simulations of channel blockers were aimed to test, as a proof of con- 
cept, a qualitative reproduction of the main experimentally AP modifications observed 
in preliminar experiments. In the case of I K „ a reduction in the maximal conductance 
simulated the application of E4031. This operation altered AP shapes quite similarly to 
what observed in the experiment, despite the fact that simulated AP prolongation was 
far less pronounced at the Early stage while at the Late stage MDP was more negative 
with respect to the experimental trace. Importantly, our model replicates the effect of 
E4031 on the frequency of spontaneous beating, which decreases in the Early stage and 
markedly increases in the Late one. 

Our experiments using 2 mM nickel were simulated by a relevant blockade of lNaCa> 
I Ca T and I Kr in the model, in accordance with the lack of selectivity of nickel at this con- 
centration and the high levels of block reported for these currents. At the Early stage 
simulations and experimental recording similarly show a residual electrical activity 
slightly different in frequency and amplitude (see Figure 6C); nonetheless, both results 
suggest that residual currents are sufficient to drive a repetitive potential oscillation of 
small amplitude occurring at depolarized potentials. For the Late stage experimental 
traces and simulations show a complete block of spontaneous activity (see Figure 6D) 
occurring at negative potential values, opening the hypothesis that in this phase, where 
other depolarizing currents such as If are downregulated, iNaCa plays a key role in sus- 
taining the spontaneous activity. Our results are consistent with those of similar experi- 
ments [45] performed with 5 mM nickel, a concentration expected to exert a higher 
level of channel blockade. In fact, with this higher nickel concentration no Ca 2+ transi- 
ents due to spontaneous activity were recorded at any developmental stage. 

Limitations 

The main limitations of our work are related to i) the shortage of data and ii) the vari- 
ability of phenotypes. We built an AP model for human embryonic stem cells derived 
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cardiomyocytes. However, many data to build the model are taken from published ex- 
perimental work on different animal species. The crossbreeding process of data from 
human and animal CMs represented a compulsory choice for our model, since some 
current properties expressed in the adult ventricular cells were unknown or not suffi- 
ciently characterized in hESC-CMs. A more accurate numerical model of "human" em- 
bryonic stem cells derived cardiomyocytes would require further experimental 
investigations. This should be considered as a first step in the effort of mathematically 
describing the ESC-CM electrophysiology aimed at capturing the qualitative 
mechanisms. 

The choice of a specific human adult ventricular AP model as the basis for the 
hESC-CM model could be perceived as a limitation of the present work. We chose as 
the "starting point" for the development of our model a modified version of the Ten- 
Tusscher model [15], in which I K „ I Ks and I Ca L were changed with respect to their ori- 
ginal formulation. From this starting point, in order to reproduce developmental 
changes in each current we modulated (through the Ral xx ratios) the current maximal 
conductances. But at the same time the kinetic of several currents were changed based 
on data from ESC-CMs: I Na , I Ca Li Ito> Iki> ^NaCa- Moreover, two relevant membrane cur- 
rents that were not present in the adult ventricular model (IcaT and If) were incorpo- 
rated in our model. All these modifications reduced the actual impact on the final 
results of the initial choice of the starting point. A sensitivity analysis with respect to 
the choice of the parent model was far beyond the scope of the present work. It is also 
worth noting that the modification of the Ten Tusscher model allowed the correct pre- 
diction of APD shortening at higher [Ca 2+ ] Q [15]. All the more recent models (e.g. [49]) 
show unrealistic APD prolongation upon increase of extracellular calcium. This aspect 
is particularly relevant for our study, in which we had to reproduce intracellular AP 
recordings obtained with 2.7 mM Ca 2+ in the external solution [50]. 

The difference in MDP between experiments and simulations is a limitation, possibly 
related to an incomplete description of the complex interconnection within the EBs. 
For example, in simulating AP recorded from EBs, only CM to fibroblast coupling has 
been considered. Other variables, such as CM to CM gap junctions likely contribute to 
the AP properties of EBs. However, such a term would dramatically complicate the sys- 
tem of differential equations and goes beyond the scope of the present work. 

Conclusions 

In conclusion, this study provides the first modelling tool able to simulate the mem- 
brane AP, the associated intracellular Ca 2+ handling properties and the modification 
occurring over the maturation process of hESC-CMs. The simulation of the transi- 
tion from Early to Late developmental stage involved: increasing I to density, declin- 
ing If contribution, decreasing iNaCa current density, Iic r decrease and In s increase. 
Moreover, an increasing contribution to Ca 2+ cycling of the Ca 2+ release from SR 
was pointed out. 

We expect to overcome inherent limitations present in the model by further ex- 
perimental investigations exploring unknown properties of basic physiology of hESC- 
CMs, possibly including different stem cell lines. Also, the combined use of novel 
pharmacological/simulated challenges will be useful to implement and validate the 
predictive potential of the model. Finally, novel challenges come from studies 
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(including drug testing) in CMs from induced pluripotent stem cells carrying genetic 
mutations [51,52]; the development of disease-specific cell lines for genetic cardiac dis- 
orders prompts toward the refinement of this mathematical modelling to address future 
directions in this field. 



Additional file 



Additional file 1: Supplementary Methods: Mathematical Modelling of the Action Potential of Human 
Embryonic Stem Cell derived Cardiomyocytes. Additional figure SI: l Ks tail current. Additional figure S2: APD 

rate-dependence [53,54]. 



Abbreviations 

AP: Action potential; APA: Action potential amplitude; APD: Action potential duration; APD X x: Action potential duration 

at XX% of repolarization; C m : Membrance capacitance; CM: Cardiomyocyte; DDR: Diastolic depolarization rate; 

EB: Embryoid body; F: Frequency of spontaneous beating; G xx : l xx current conductance; hESC: Human embryonic stem 

cell; hESC-CM: Human embryonic stem cell-derived cardiomyocyte; l bCa : Background Ca 2+ current; l CaL : L-type Ca 2+ 

current; l CaT : T-type Ca 2+ current; l f : Hyperpolarization activated funny current; l Kr : Rapid delayed rectifying K + current; 

1^: Slow delayed rectifying K + current; l K1 : Inward rectifying K + current; l| eak : Leakage Ca 2+ current; l maxX x: U maximal 

current; l Na : Na + current; l NaCa : Na + /Ca 2+ exchanger current; l NaK : Na + /K + pump; l pCa : Sarcolemmal Ca 2+ pump; 

l re i: Release Ca 2+ current; l to : Transient outward K + current; l up : Uptake Ca 2+ current; MDP: Maximum diastolic potential; 

Ra !xx : Variable fraction of the l xx current maximal conductance in the adult model; SR: Sarcoplasmic reticulum; 

Vmax: Maximal upstroke velocity of the membrane potential; VmaxCa.upstroke: Maximal upstroke velocity of the Ca 2+ 

transient; VmaxCa.decay: Maximal decay velocity of the Ca 2+ transient. 
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